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1. INTRODUCTION 



There are several recently proposed classes of 

empirical probability density function [1,4, 5, 7] all 

generally considered to be superior to the classical 

histogram estimates. The class considered in this paper is 

based on independent observations, i.e. X ,X ,...,X are 

12 n 

independent and identically distributed random variables 

with continuous unknown density function f (x) . The method 

used to estimate f (x) is that proposed by Rosenblatt; 

denoting the estimate by f (x) , we define 

n 



f (x) 
n 



1 

nBTnr 




r x - X 
WI 



L "57n F 



a] 



r 



where W (u) is a bounded non-negative integrable weight 
function with 




W (u) du 



1, 



and b (n) is a positive bandwidth function which tends to 
zero as n --> °° , but is such that o[ b (n) ] = 1 / n. Thus we 



might have b (n) 'w 



- 1/2 

n , for example. 



We note that 
density functions 



all estimates of this form are 
for a given set of observations 



themselves 
that is, 



f (x) > 0, 
n 




(x) dx 



1 . 



Since the 
parameter 
non-s ta tiona 



X ’s are random variables, f 
j n 

stochastic process, but 
ry. 



(x) 

it 



is a 
is 



continuous 

clearly 



1 



The estimate f (x) can be shown to be locally biased 
n 

for any value of x under relatively mild conditions [4]. 

Our object in this paper is to investigate a global measure 

of how good f (x) is as an estimate of f (x) . The measure 

n 

was originally proposed by Bickel and Rosenblatt [2] and is 



given by 



6 <n> - f 



j-f (X) - f (X) -J 2 

— - n*r 



dx . 



Since the value of £ (n) will vary with each realization of 

X ,...,X , it is a statistic or function of the n random 

1 n 

variables. A possible application for such a statistic 

would be in good ness -of -fit type tests, in an analogous 
manner to the more familiar Kolmogorov-Smirnov test. 



Bickel and Rosenblatt [2] have established that if 
-2/9 

b(n) = o[ n ] as n --> 00 and if a (x) is a bounded. 



piecewise smooth integrable function then 

b (n) ^ [ nb (n) J " [f (x)-f(x)] 2 a(x)dx - J " E (x) a (x) dx j" ( z) 2 dz ] 

is asymptotically normally distributed with zero mean and 
variance 



2W<*>(0) J* a(x) 2 f(x) 2 dx , 



as n --> 00 , where W <A >(0) is the fourth convolution of K 

with itself. Thus, B (n) has an asymptotically normal 
distribution, regardless of the underlying density f (x) . 



A problem in this situation is that, unlike the 
Kolmogorov-Smirnov test statistic, the statistic B (n) is 



2 



not distribution-free. Further, its exact distribution for 
any finite value of n does not seem to be mathematically 
tractable. We thus examined some representative cases 
through simulation, hoping that 3 (n) would be fairlv robust 
with rapid convergence to the asymptotic distribution. It 
was also hoped that the simulations would cast light on 
these conjectures and perhaps suggest some unexpected 
results. 
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2. SIMULATION 



The primary object of the simulation was to investigate 
the distribution of the statistic 8(n) : 

r c f (X) - f (X),2 

6 (n) = / L n i dx f 

over a suitable range of integration. We performed 
simulations with synthetic sampling from both uniform and 
Cauchy distributions; the triangular weight function 

/ 1 ~ I u | , if I u | < 1 

W(u) = < 

l 0 , otherwise 

was used to evaluate f (x) in both cases. We found little 

n 

difference as far as 8(n) was concerned between the 

triangular and other "smoother” (e.g., quadratic) weight 
functions for our samples of from 100 to 1500 deviates. 

A. UNIFORM RANDOM VARIABLES 



In the case of uniform (0,1) random variables, we have 



f (x) 




0 < x < 1 
otherwise . 



Thus, g (n) becomes. 



3 (n) 



= /I- b(n) 

J b (n) 



[f (x) 

n 



1]2dx . 



( 2 . 1 ) 



The limits of integration are from b(n) to 1 - b ( n ) 

instead of from 0 to 1 to avoid the marked bias of f (x) 

n 

near 0 and 1. As long as b (n) < x < 1-b(n) , though, f (x) 

n 
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is unbiased: 



E[f (x) 1 
n 



= 1 / W r x-y , dy 

B7nyy 0 *-57np 



EW. 



/*x + b (n) 

/ r 1 - 1 X "ZJ -1 dy 

nryx-b(n) L 5TnP 



Elnf 



r/' x+b(n) dy - f x 

nf L J x-b(n) J x- 



b (n) 



x- y dy 

Elnr 



x + b (n| 



■/ 



y-x dy 1 

B15T 



1 r 2b (n) - 1 LkinLJ£ - _1 LklnlJi 

BlnB L B^nf 2 B^nf 2 J 



= 1 



Also, for the same range of x, 

r p. r x “ X t t 

Var[f (x)] = Var J 1 J W| 3| J 

n t nB^nT 3=1 L oT n T J J 



1 

n^BtnJ* 7 j 



^ Var w[ X jl 
= 1 



rX - X -« 
= 1 Var W i 

nBTnF 7 <-"5Tnr J 



= 1 \ T H 2 rX-y; t dy - \ f R r x-y i dy ] ] . 

nETnF 7 Fy o ^57nr 1 U 0 «-5 Tnf J -» J 



0 L B7np 



0 L E7np 



Since f (x) is a piecewise linear function when a 
n 

triangular weight function is used, the integral in (2.1) 

can be evaluated in principle but the work becomes 

prohibitive for even moderate sample sizes. We thus 

approximated the integral using Simpson's rule with 100 

equal subintervals. The results were found to be 

satisfactory in the sense that the value did not change 

appreciably when a finer grid (up to 500 subintervals) was 

used. In general, we found that a larger sample size 

required a finer grid; apparently the value of f (x) changes 

• n 
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more rapidly over a small interval when n is large. 



3 



/ 



We used 

1/2 

n , 



three different 
1/2 

1 / n and 1 



bandwidths in the uniform case: 
/ n. For each bandwidth sample 



sizes of 100, 200, 500, 
that a total of 15 
experiment consisted of 
of which resulted in 
g (n) using (2.1) . The 
were divided into five 
that variability of the 
between sections. 



1000 and 1500 were investigated so 
experiments were carried out. Each 
2000 independent replications each 
the calculation of a single value of 
replications for a given experiment 
sections of 400 observations each so 
simulation results could be assessed 



Besides the 400 observed values of g (n) , the computer 
output for each of the 75 sections included a histogram, an 
empirical log-survivor function plot, an empirical CDF plot 
and a normal probability plot. A histogram and an empirical 
log-survivor plot were also computed for the pooled sample 
of 2000 for each experiment. These plots are all reproduced 
in reference [3]; some of the more interesting cases are 
included in Section 4. 



It was found that a better picture of the distribution 
of the data resulted when the empirical density function of 
the g (n) 's was plotted over the histogram plot. A fairly 
wide bandwidth was needed to suppress large fluctuations in 

1/2 

f (x) ; it was found that b(n) = R / n was a fairly robust 



n 

choice. (R denotes 
minimum value] of the 
Figures in Section 4 
this bandwidth and the 



the sample range [maximum value - 
B (n) sample.) The solid lines in the 
are empirical density estimates using 
triangular weight function. 
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B 



CAUCHY RANDOM VARIABLES 



The Cauchy density function is 



f(jc) = 1 

!rTT+x*r 

We used the same density estimator as in the uniform case: 

n |-X * X i 

f (X) = 1 ^ H| j| , 

n nBTnf j = 1 <-“5Thr - 1 

and again the triangular weight function. We chose a range 
of integration (-3, + 3): 



6 



f+3 r f (x) - f (x) p 

(n) = I n 

J -3 n*r 



dx 



This range comprises 80S of the probability mass for this 
distribution. Again, Simpson's rule was used to approximate 
the integral; in this case a grid of 600 subintervals was 
selected after examining 100, 300, 600 and 900 subinterval 
grids . 



The Cauchy distribution was chosen because for finite n 

f (x) has a bias component; this component usually decreases 
n 

with bandwidth for a fixed value of n, although the 

pointwise variance of f (x) increases with decreasing 

n 

bandwidth. It seems likely that the variance of g (n) would 



also decrease under these conditions, as indeed it was 
observed to do. 



Three bandwidths were also employed in the Cauchy case: 
1/2 1/2 1/2 

1 / n , 3 / n and 20 / n , the last one representing 

a case in which bias in the estimator f (x) plays a major 

n 
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role in 



the distribution of 



6 (n) . 



The same five 



sample 



sizes were used here for each 
uniform simulations; output 
experiments was obtained just 



bandwidth as were used 
from the fifteen 
as in the uniform case 



for the 
Cauchy 
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3. TABULAR RESULTS AND GAMMA FITS 



Using the asymptotic result obtained by Bickel and 
Rosenblatt [5], for a uniform random variable the quantity 



- 1/2 

b(n) {nb(n) 



I 



1-b (n) 
b (n) 



|f (x) - 1 1 2 dx 
n 



1 - 2b (n) 



/ 



W (U) 2du} 



is asymptotically normally distributed with mean 0 and 
variance 



2W<*> (0) [ 1-2b(n) ] 



as n — > <» if nb{n) — > <*> and b (n) 



-2/9 

o (n ) . 



For the 



triangular weight 



function , 




(u) 2 du 



2 

3 



and W< 4 >(0), the fourth convolution 
zero, is 302/630. 



of W with itself at 



Var[ 



From the above 
Er 6 (n) ] = z[ 

g(n) ] = Var 




expres 
f 1-b(n 

J b(n) 



1-b (n) 



b(n) 



If 



sions. 

If (X) 
n 

(x) -1 
n 



we get 




1 J 2 dx 1 

J 


^ 2 1-2b(n) 

3 _ nE7n) 


2 dx ] -v 


2HC4)^M-2binU 



Comparisons of the simulated values for the uniform 
experiments with the conjectured ones are tabulated in Table 
III.1 (means) and Table III. 2 (variances). Especially for 
small bandwidth the agreement between the asymptotic and 
simulated variances is very good even for small n (n = 100) . 
The same is true for expected value, although convergence is 
slower than for the variance and again slower for large 
bandwidth. 
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table 1 1 1 - 1 Comparison of estimated mean values and asymptotic 
mean values of B(n) for different bandwidths and sample sizes. 



n 


b(n) = 3//n 


E (B(n) ) 


E (S (nil / ( 1 - 2b (si) J 


Conjectured 


Computer output 


100 


. 3000 


.0089 


.0222 


.0127 


200 


.2121 


.0090 


.0157 


.0109 


500 


.1342 


.0073 


.0099 


.0075 


1000 


.0949 


.0057 


.0070 


.0058 


1500 


.0775 


.0048 


.0057 


.0051 




b (n) = l/»^n 








100 


.1000 


.0533 


.0667 


.0583 


200 


.0707 


.0405 


.0471 


.0415 


500 


.0447 


.0271 


.0298 


.0269 


1000 


.0316 


.0197 


.0211 


.0197 


1500 


.0258 

i 


.0163 


.0172 

i 


.0168 

i 



TABLE 1 1 1. 2 Comparison of estimated standard deviation values and 
asymptotic standard deviation values of 8(n) for different band- 
widths and sample sizes. 



n 


b (n) =3/ /rT 


a(B(n)) 


a (B (n) ) / ( 1 -2b (n) ) 


Conj ccturc J 


Computer output 


100 


. 3000 


.0113 


.0283 


.0115 


200 


.2121 


.0081 


. 0141 


.0088 


500 


. 1342 


.0046 


.0063 


.0047 


1000 


.0949 


.0029 


. 0036 


. 0030 


1500 


.0775 


.0022 


. 0026 


.0023 




b(n) = l//n 








100 


. 1000 


.0277 


.0346 


.0315 


200 


.0707 


.0171 


.0199 


.01S9 


500 


.0447 


.0088 


.0097 


.0092 


1000 


.0316 


.0053 


.0057 


.0056 


1500 


.0258 . 


.0040 


.0042 


.0043 



lO 



In contrast to the moments, the distribution of & (n) 
converges very slowly. The complete results (reference 
[3]) reveal that the histograms and empirical density 
functions of the 8(n)'s are all skewed to the right; see 
Figures IV. 1 to IV. 9 for examples. 

The form of the histograms as well as the log-survivor 
plots suggested that the B(n) statistic is approximately 
Gamma (9, k) distributed, where the Gamma density is given by 



k-1 -x/9 




and the mean and variance are 

E[ X ] = k9; 

Var[ X ] = k9 2 . 

Accordingly, estimates ?. and 9 of k and 9 for each 

experiment were obtained from the sample of 2030 6 (n) 's. 

Shenton and Bowman's almost unbiased estimators for the 
Gamma distribution [6] were used; these give reasonable 
results when k > 0.5, as in this case. The estimate values 
are tabulated in Taole III. 3; also tabulated are estimates 

of the standard deviation of K and 9 which were obtained 

from the five sections in each experiment. A parametric 
density estimate is thus obtained for the £ (n) sample; it 

may be compared with the non-parametric estimate f (x) by 

n 

examining the graphs in Section U, where the Gamma density 



function is plotted with a dashed line. 



for 



Fitted 



Gamma 



TABLE III. 3 
Distribution 

DISTRIBUTION 



UNIFORM 



CAUCHY 



Estimated 
»r 6 (n) . 


Parameters 


b(n) 


n 






100 


3 




200 


5 


1 / /n 


500 


8 




1000 


13 




1500 


17 




100 


1 




200 


1 


3 / /n 


500 


2 




1000 


4 




1500 


5 




100 


40 




200 


39 


1 / n 


500 


33 




1000 


32 




1500 


31 




100 


22 




200 


3 2 


1 / /n 


500 


60 




1000 


79 




1500 


101 




100 


9 




200 


12 


3 / /n 


500 


20 




1000 


29 




1500 


34 




100 


7 




200 


4 


20 / /n 


500 


3: 




1000 


4' 




1500 


5; 



± 0.206 
. 780 
± 0.659 
. 881 
± 0.899 
. 011 
± 0.796 
.316 
± 1.467 

. 153 
± 0.048 
.718 
± 0.174 
.707 
± 0.241 
. 028 
± 0.241 
. 248 
± 0.423 

. 337 
± 2.555 
.511 
± 2.347 
. 649 
± 1.820 
. 033 
± 3.305 
.712 
± 1.999 

. 362 
± 1.488 
. 305 
± 2.022 
. 147 
± 4.009 
. 897 
± 6.608 
. 100 
± 7.783 

.272 
± 0.406 
. 744 
± 0.645 
. 701 
± 1.673 
. 303 
± 2.541 
. 265 
± 3.963 

. 103 
± 0.217 
. 144 
± 0.069 
.445 
± 0.161 
.211 
± 0.357 
. 385 
± 0.335 



B 



0 .01390 
± .00095 
0 . 007 15 
± .00095 
0.00311 
± .00029 
0.00153 
± .00009 
0.00095 
± .00008 

0 . 00967 
± .00058 
0.00588 
± .00078 
0.00281 
± 0.00026 
0.00145 
± .00007 
0.00096 
± .00008 

0 .01616 
± .00117 
0.01675 
± .00106 
0. 01953 
± .00111 
0. 02059 
± .00244 

0.02089 
± .00124 

0.01745 
± .00114 
0.00864 
± .00054 
0 . 00293 
± .00022 
0. 00157 
± .00014 
0.00102 
± .00007 

0.01331 
± .00062 
0.00709 
± .00037 
0 . 00277 
± .00022 
0.00140 
± .00012 
0.00099 
± .00010 

0. 00776 
± .00035 
0 . 00619 
± .00009 

0. 00312 
± .00016 
0.00152 
± .00009 
0. 00095 
± .00005 
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4. GRAPHICAL RESULTS AND GENERAL DISCUSSION 



The graphs for the following experiments have been 
reproduced from [3] because they give the greatest insight 
into the distribution of g (n) ; these graphical results are 



more informative than the tabulated 
Gamma fits of the previous Section. 



means, variances and 



Figure 


Random Variable 


n 


b(n) 

1/2 


Jc 


4. 1 


Uniform 


200 


3/n 

1/2 


1.718 


CN 

• 

^3* 


Uniform 


500 


1/n 

1/2 


8.881 


4.3 


Uniform 


1 500 


1/n 


17.316 


4.4 


Uniform 


200 


1/n 

1/2 

1/n 

1/2 

3/n 

1/2 


39.511 


4.5 


Cauchy 


100 


22.362 


4.6 


Cauchy 


100 


9.272 


4.7 


Cauchy 


1500 


2 0/n 

1/2 

3/n 

1/2 


5. 385 


00 

• 


U nif orm 


1500 


5.248 


4.9 


Uniform 


100 


1/n 


3.969 


In interpreting the graphs 


we can be 


guided 


by c 


ristics. 


In the case of 


a density 


estimate 


f (X) 



bandwidth b(n) there is dependence within a range 



n 

of 



order 



b(n) and an approach to independence for points separated by 
a distance of order larger than b (n) . Thus in the case of 
uniform random variables the integral g (n) could be thought 
of as having the equivalent of the order of [1-2b(n) ]/b(n) 
independent summands. In the first case (Figure 4.1; n=200, 

b(n)=3//n, Tt=1.7l3) we obtain 



(1 - 3/2/10) / [3/(1 0/2 ) ] = 2.71 . 

This is rather small so that one does not expect a good 
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Gaussian fit. We give X from the previous Section since 2Jc 

may be interpreted as an equivalent number of degrees of 

freedom; the larger the fitted k, the closer ne are to 

normality. In a loose sense it is clear that a gamma fit is 
likely to be more appropriate and this is confirmed by 
looking at the graphs. 

In the second case (Figure 4.2; n=500, b(n)=1//n, 

k=8.881) ve have 

(1 - /2/10) 10/2 = 12.14 , 

which is a bit larger. It is interesting to note that the 
estimated (smoothed) density function of e(n) gives us 
greater insight apparently in all cases. Here we see the 
beginning of an approach to asymptotic normality though it 
is still suggested that a Gamma fit might be appropriate. 

The next case (Figure 4.3; n=1500, b(n)=1//n, k=17.316) with 

[ 1 - 1/ (5/T5) ] 10/To = 36.73 

shows a closer approach to normality. It may be seen that 

the major departure between the parametric and 

non-parametric density estimates occurs in the vicinity of 

the mode where f (x) tends to fluctuate about the true 

n 

value. The fit in the tails appears excellent in all cases. 



The 
k=39. 511) 
suggested 



next uniform case (Figure 4.4; n=200, b(n)=1/n, 
is strictly speaking outside the range of results 
by the paper of Bickel and Rosenblatt [2]. Here 



f (x) is asymptotically compound Poisson rather than 
n 
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asymptotically normal. Nonetheless we notice that it looks 
as if a Gaussian fit would be very good and this is 
consistent with the magnitude of our crude index 

(1 - . 01) 200 = 198 . 

It would be interesting for someone to prove the suggested 
asymptotic normality. 

In the simulation of sampling from a uniform 

distribution, the density estimator has no bias. To 
investigate the effect of bias, we repeated the uniform 
experiments for Cauchy-distributed random variables, 
integrated over the range -3 + b(n) to 3-b(n). The first case 

(Figure 4.5; n=100, b(n)=1//n, K=22. 362) has index 

(6 - .2) 10 = 58 

and one notices that a Gaussian fit looks very good. The 

next case (Figure 4.6; n=100, b(n)=3//7T, Jt = 9.272) has index 

(6 - .6) 10/3 = 17.66 , 

and a Gaussian fit looks fair but not good. In the last 
Cauchy case one expects substantial bias (Figure 4.7; 

n=1500, b(n)=20//n, K=5.385) and the crude index is 

(6 - 4//T5) ST5/2 = 9.6 2 . 

A Gamma fit is suggested. Altogether the effects of bias 
don't seem to be that extreme when sampling from the Cauchy 
distribution but this may be due to the fact that the Cauchy 
density is a very smooth function. 

The last two cases involve sampling from the uniform 
distribution again but with different sample sizes and 

bandwidths. Figure 4.8 is for n=1500, b(n)=3//n and 

K=5.248, while Figure 4.9 is for n=100 and b(n)=1//n for 

which Tc=3.969. 
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The problem in using 8 (n) as a measure of gaodness of 
fit in the non-limiting 3arama case is to determine k and 9. 
If one wishes to fit the Samma distribution using the method 
of moments, one can use the fact that the mean and variance 
of 8 (n) should be approximately (on asymptotic grounds) 
W< 2 >(0) and b(n)W<*>(0), respectively. One might then use 



as estimates of k and 9. The results in Section 3 suggest 
that this procedure should produce adequate results except 
when there is appreciable bias in the density function 
estimate . 




0 
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Figure 4.2. Distribution of the statistic B 

uniform random variable with n = 500 and bandwidth 

The solid line shows the Rosenblatt empirical 
function of the 8 (n) 's while the dashed line is 

Gamma density function with Ic = 8.881 and 9 = . 003 1 
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Figure 4.3, 



Distribution of the statistic 3 (n) for 



uniform random variable with n = 1500 and bandwidth 1 / y/H. 
The solid line shows the Rosenblatt empirical density 
function of the g (n) ’ s while the dashed line is a fitted 



Gamma density function with 1c = 17.316 and 9 = .00095. 
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Figure 4.4. Distribution of the statistic 8(n) for a 
uniform random variable with n = 200 and bandwidth 1 / n. 
The solid line shows the Rosenblatt empirical density 
function of the 8 (n) 's while the dashed line is a fitted 

Gamma density function with 1c = 39.511 and 0 = .01675. 
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Figure 4.5. Distribution of the statistic 6 (n) for a 

Cauchy random variable with n = 100 and bandwidth 1 / x/n. 

The solid line shows the Rosenblatt empirical density 
function of the 3 (n) 's while the dashed line is a fitted 

Gamma density function with Tc = 22.352 and 9 = .01745. 
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Figure 4.6. Distribution 

Cauchy random variable with n 

The solid line shows the 
function of the 6(n) 's while 

Gamma density function with Tc 



of the statistic 6 (n) for a 

= 100 and bandwidth 3 / \/n. 

Rosenblatt empirical density 
the dashed line is a fitted 

= 9.272 and 0 = .01331. 
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Figure 4.7. Distribution of the statistic g (n) for a 

Cauchy random variable with n = 1500 and bandwidth 20 / ^/n. 
The solid line shows the Rosenblatt empirical density 
function of the $ (n)*s while the dashed line is a fitted 

Gamma density function with K = 5.385 and & = .00095. 
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Figure 4.8. Distribution of the statistic 3 (n) for a 

uniform random variable with n = 1500 and bandwidth 3 / y'n. 
The solid line shows the Rosenblatt empirical density 
function of the g (n) 's while the dashed line is a fitted 

Gamma density function with 1c = 5.248 and 0 = .00096. 
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Figure 4.9. Distribution of the statistic B (n) for a 

uniform random variable with n = 100 and bandwidth 1 / v/n. 

The solid line shows the Rosenblatt empirical density 
function of the B (n) 's while the dashed line is a fitted 

Gamma density function with 1c = 3.969 and 0 = .01390. 
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